1/1 


AD-A17S  190 


UNCLASSIFIED 


RATIONAL  ARITHHETIC  IN  FLOATING-POINT<U>  CALIFORNIA 
UNIV  BERKELEY  CENTER  FOR  PURE  AND  APPLIED  HATHEHATICS 
H  KAHAN  SEP  86  PAN-342  N00014-8S-K-0180 


A175  190 


CENTER  FOR  PURE  AND  APPLIED  MATHEMATICS 
UNIVERSITY  OF  CALIFORNIA,  BERKELEY 


PAM-343 


RATIONAL  ARITHMETIC  IN  FLOATING-POINT 


W.  Kahan 


Q 

< 


September  1986 


This  report  was  done  with  support  from  the  Center  for 
Pure  and  Applied  Mathematics.  Any  conclusions  or 
opinions  expressed  in  this  report  represent  solely 
those  of  the  author(s)  and  not  necessarily  those  of 
the  Center  for  Pure  and  Applied  Mathematics  or  the 
Department  of  Mathematics. 


DISTRIBUTION  STATEMENT  A 

Approrod  for  public  relea»«; 
Distribution  Unlimited 


RATIONAL  ARITHMETIC 


FLOATING-POINT 


$ 


ELECTE 
DEC  1  6  1986 


W.  Kahan 

Math.  Dept.,  and  E.  E.  St  Computer  Science  Dept, 
University  of  California  at  Berkeley 
Sept.  20,  1986 


Abstract:  Calculating  M/N  1*  A/B  £  C/D  in  lowest  terms,  given 

the  integers  A,  B ,  C  and  D,  is  a  task  taught  in  Elementary 
schools;  and  it  is  an  easy  exercise  in  Computer  Programming 
too  provided  the  given  integers  must  be  less  than  half  as  wide 
as  the  widest  integers  that  can  be  handled  conveniently  by  the 
computer's  hardware  or  by  its  programming  language.  But  that 
program  becomes  much  more  complicated  (and  slower)  if  it  is 
naively  expected  to  perform  correctly  whenever  all  six  of  our 
integers  A ,  B ,  C,  D,  II  and  N  are  allowed  to  grow  almost  as 
wide  as  those  widest  convenient  integers.  This  simple  task 
illustrates  why  the  art  of  programming  entails  sometimes  a 
delicate  balance  between,  on  the  one  hand,  the  simplicity 
and  aesthetic  appeal  of  the  specif ications  and,  on  the  other 
hand,  the  complexity  and  efficiency  of  the  implementation. 

r 

Introduction: 

The  obvious  way  to  calculate 

M/N  1=  A/B  ±  C/D  in  lowest  terms 
is  to  first  calculate 

M*k  :=»  A*D  ±  B*C  and  N*k  :  -  B*D 
and  then  divide  them  by  their  Greatest  Common  Divisor 

k  .*=  gcd  (M*k,  N*k)  . 

But  the  obvious  way  is  no  way  to  calculate 
31  /  1897  51872  =  1234  56799  /  123456  -  9882  97396  /  988291 

on  a  calculator  that  carries  only  ten  significant  decimals 
because  first 

M*k  :=  1234  56799  *  988291  -  123456  *  9882  97396 

*  12201  12433  40509  -  12201  12433  20576 

=  1 9933 

and 

N*k  :=  123456  *  988291  *  12  20104  53696 

would  have  to  be  calculated  in  order  to  reveal 

k  1=  gcd  ('  1 9933 ,  12  20104  53696.)  =  643  . 

On  that  calculator,  the  two  fifteen-digit  products  would  both 
round  to  the  same  value  (  12201  12433  00000  )  to  ten  significant 
digits,  yielding  zero  for  M*k  ;  and  N*k  would  get  rounded  off 
too.  However,  because  the  desired  final  results  M  =  31  and 
N  =  189751872  can  be  held  exactly  in  that  calculator,  a  wav  to 
compute  them  exactly  ought  to  exist.  An  algorithm  that  does  so 
without  merely  simulating  arithmetic  to  at  least  fifteen  Digits  is 
the  subject  of  this  note.  The  algorithm  is  not  simple,  but  it  is 
far  simpler  than  simulating  multi-word  arithmetic  in  BASIC  . 


The  Computing  Environment: 

There  are  limits  to  the  widths  of  the  inteqers  and  floating-point 
variables  supported  conveniently  in  programming  languages  like 
Fortran,  BASIC,  Pascal  and  C.  Integers  on  some  computers  may 
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be  no  wider  than  16  bits,  running  iron  -32768  to  32767;  on 
most  other  computers  the  integers  occupy  32  bits,  running  -from 
-21474  83648  to  21474  83647  .  Integers  bigger  than  that  lose 
their  leftmost  bits  to  Overflow,  usually  without  any  warning 
accessible  to  the  higher — level  language  program.  Floating-point 
variables,  limited  to  24  significant  bits  on  some  machines,  to 
53  or  36  on  most  others,  can  handle  much  bigger  integers;  but 
integers  bigger  than 

2.0 **  -  167  77216.0  or 
2.0”  =  9  00719  92547  40992.0  or 
2.0”  =  72  05759  40379  27936.0  respectively 
lose  their  rightmost  bits  to  Roundoff ,  and  consequently  become 
multiples  of  powers  of  2  even  when  ideally  they  should  have  been 
odd.  Similarly,  on  a  typical  ten-digit  calculator,  integers 
bigger  than  1  00000  00000  get  rounded  off  to  multiples  of  powers 
of  ten.  Rounding  errors  occur  without  any  warning  to  the  program 
(except  on  machines  that  conform  to  IEEE  standards  754  and  p854, 
which  require  that  rounding  errors  signal  Inexact.)  That  lack  of 
warning  obliges  programmers  to  clutter  some  programs  with  tests  of 
the  magnitudes  of  all  intermediate  results  lest  incorrect  final 
results  be  produced  with  no  indication  that  they  are  wrong. 

Let  A  stand  for  the  smallest  positive  integer  beyond  which  some 
digit  must  be  lost  to  overflow  or  roundoff;  the  previous 
paragraph  tenders  values  of  A  appropriate  for  various  machines. 

A  is  what  is  meant  by  "the  widest  integer  that  can  be  handled 
conveniently  by  the  computer's  hardware  or  by  its  programming 
language."  The  obvious  way  to  calculate  M/H  described  above 
would  obviously  work  if  |4*0|,  |S*C{  and  10*01  were  all  somewhat 

smaller  than  A  ,  as  would  surely  be  the  case  if  Ml,  |S|,  | C I 
and  | D I  were  all  somewhat  smaller  than  /A  .  The  vagueness  here 
implied  by  the  word  "somewhat"  allows  for  sloppy  implementations 
of  floating-point  arithmetic  that,  on  some  machines,  introduce 
unnecessary  rounding  errors  when  integer  results  approach  A  too 
closely.  Notwi thstandi ng  that  vagueness,  an  algorithm  will  be 
presented  that  calculates  M  and  N  exactly  whenever  they  and  the 
given  integers  A,  B ,  C  and  D  are  all  somewhat  smaller  in 
magnitude  than  A  rather  than  merely  /A  . 


Rem,  gcd,  and  Lowest  Terms: 

Our  algorithm  will  require  certain  utilities  which,  if  not 
already  present  in  the  programming  environment,  will  have  to  be 
programmed  from  scratch.  Reducing  (M*k) / (N*k)  to  its  lowest 
terms  M/N  requires  that  k  *  gcd (M*k,  N*k.)  be  computed;  and 
the  fastest  ways  to  compute  gcd's  require  that  remainders  be 
computed.  Let 

remM,  y)  .*  =  x  -  y*( the  integer  nearest  x/y)  provided  y  f  0  .  4* 

This  is  consistent  with  the  definition  of  the  operation  rern  that 
must  be  present  ir  programming  environments  that  conform  to  the 
IEEE  standards  754  and  p834  for  floating-point  arithmetic.  In 
other  environments,  rem  must  be  composed  from  other  primitives. 

In  Fortran  the  generic  intrinsic  function  MOD  (for  INTEGER'S,  - 
AM0D  for  REALs,  DM0D  for  DOUBLE  PRECISION)  serves  to  define  REM 
thus:  .  . 

/ ».  m  keri&- 
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GENERIC  FUNCTION  REMlx,  y) 

REM  •  MOD <  x ,  y) 

IF  (  ABS (REM)  .GT.  ABS(y-tf£7f)  )  REM  -  y  -  REM 

RETURN 

END 

Absent  rent  and  MOD,  the  -following  procedure  might  be  used: 
function  remix ,  y)\ 
q  : *  x/y  ; 

n  : *  q  rounded  to  the  nearest  integer  : 

return  rem  x  -  y*n  ; 

end. 

Both  procedures  can  malfunction  when  x  approaches  or  exceeds  A 
in  magnitude;  the  following  example  will  show  how  roundoff  in 
x/y  and  y*n  causes  trouble. 

Suppose  floating-point  arithmetic  is  rounded  to  six  significant 
decimals,  for  which  A  =  1000000  .  Now  take  x  -  999999.0  and 
y  =  9901.0  ,  whereupon  x/y  53  100.99979  80002  ...  must  round 
to  q  =  101.000  .  Then  n  »  <7  ,  but  y*n  -  A+l  must  round  to 
A  ,  which  wrongly  returns  -1.0  instead  of  -2.0  for  rem  . 
Similar  rounding  errors  inside  the  implementation  of  AM0D  can 
return  -1.0  instead  of  9899.0  for  AM0D (999999. 0 ,  9901.00)  . 

If  the  quotient  x/y  were  chopped  instead  of  rounded,  no  such 
malfunctions  could  occur.  With  rounding,  they  can  be  avoided  bv 
keeping  I  .v  I  and  |y|  both  smaller  than  A/2  .  If  the  error  bound 
for  floating-point  division  is  vague,  as  it  is  for  CRAYs .  we 
can  compensate  for  ignorance  by  further  restricting  | x I  and  |y|  : 
that  is  why  phrases  like  "somewhat  smaller  than  A  "  have  been 
uttered  above. 

Having  found  a  way  to  compute  rem/x,  y)  well  enough  that 
(x  -  rem(x,y))/y  is  an  integer  exact l y,  and 
I  rem('.v,y.)  I  <  I  y  I /2  roughly, 

we  may  use  it  to  compute  Greatest  Common  Divisors  quicklv  thus: 
function  gcdix ,  y)t 

while  y  5*  0  do  <  temp  :=  y  ; 

y  remix,  y)  ; 
x  .‘  =  temp  }  ; 

return  gcd  :=  I x I  ; 
end. 

Besides  the  usual  properties  for  positive  integers  x  and  y  , 
namel y 

gcd ix , y )  is  the  largest  integer  such  that 
x/gcd/x,y.)  and  y/gcd('x,y\>  are  both  integers  exactly, 
this  procedure  gcd/x,  y.)  has  useful  properties  when  its 
arguments  are  negative  integers  or  zero; 

gcd  ix,y)  »  gcd /lx  I,  I  y  I  .)  and  gcd/x,0.>  =  gcd/0,x;  *  |x|  . 

These  properties  simplify  the  explanation  of  the  assertion 

"  M/N  is  in  lowest  terms,  " 

which  shall  now  be  taken  to  mean  that  integers  M  and  N  satisfy 
"  N  >  0  ,  and  either  gcd iM,N>  =1  or  M  =  N  =  0  . 

We  shall  abbreviate  "in  lowest  terms"  to  " ilt"  and  use  it  not 
only  as  an  adjective  but  also  as  an  operator  that  maps  pairs  of 
integers  to  pairs  thus: 
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■function  llt(x ,  y)i 

g  :=  copysignf  max  -Cgcd  ( x ,  y )  ,  13  ,  y  )  ; 
return  lit  ‘.=  (  x/g ,  y/g  )  ; 

end. 

Now  asserting  that  (M,  N )  *  Iltf*,  y)  means  the  same  thing  as 

M/N  -  x/y  ilt  . 


Idealized  Rational  Operations 

The  mapping  Ilt  provides  a  unique  pair  o-f  integers  (M ,  N)  to 
represent  each  rational  number  M/N  ®  x/y  ilt  ,  including  also 
+  1/0  -  +00  ,  as  well  as  a  representati on  -for  the  entity  0/0 
called  "NaN"  (-for  "Not  a  Number")  in  the  IEEE  standards  -for 
-floating-point  arithmetic.  But  those  standards  also  speci-fv  how 
+0  and  -0  will  behave  ari thmet i cal  1 y  in  case  a  programmer  chooses 
to  distinguish  them,  something  that  cannot  be  done  use-full y  on 
most  machines  that  do  not  con-form  to  those  standards.  Without  a 
well-behaved  signed  zero,  attempts  to  distinguish  between  +oo 
would  run  afoul  o-f  identities  like  M/N  -  -1  / (-N/M)  when  11  =  1 
and  N  =  0  .  That  is  why  we  shall  herein  regard  oo  as  unsigned, 
like  0  ,  as  if  the  ends  o-f  the  real  axis  had  been  lifted  and 
joined  to  form  a  circle  out  of  it.  Rational  operations  consistent 
with  that  picture  are  defined  in  a  familiar  way  as  follows: 

A/B  +  C/D  :=  ( A*D  +  B*C) / (B*D)  ilt  respectively  ; 

(A/B)  *  (C*D)  .*=  (A*C) / (B*D)  ilt  ; 

(.A/B)  -r  (C/D)  :=  (A*D)/(B*C)  ilt  ; 

A/B  -  C/D  just  when  A*D  =  B*C  but  I A*C |  +  IB#PI  #  0  . 
Thus,  the  set  of  all  rational  numbers,  augmented  by  oo  and  0/0  , 
constitutes  a  system  closed  under  the  rational  operations  so 
defined.  But  the  subset  of  rational  numbers  M/N  representable 
conveniently  on  our  computer,  those  for  which  I M I  and  IW[  do 
not  exceed  A  ,  does  not  constitute  a  closed  system:  instead  it 
poses  a  challenge  to  implement  the  rational  operations  correctly 
for  those  operands  and  results  that  do  lie  within  the  subset. 

Implementations  of  multiplication,  division  and  equality-testing 
are  entirely  strai ghtf orward ,  as  follows  below,  provided  all 
operands  are  ilt.  In  other  words,  the  operands  are  presumed  to 
be  pairs  of  integers  that  will  pass  unchanged  through  the  function 
Ilt  ,  and  the  results  will  do  the  same  provided  their  magnitudes 
are  somewhat  smaller  than  A  . 

function  Product  (A,  B ,  C,  D)  :  ...  to  get  (A/B)*(C/D'>  lit 

k  :=  max C 1 ,  gcd(A,  D) 3  ;  m  : =  max  Cl,  gcdiS,  C ) 3  ; 

return  Product  .‘ =  (  (A/k )  *  (C/m)  ,  (B/m)*(D/k'>  )  ; 

end. 

Note  that  if  the  final  results  are  all  somewhat  smaller  in 
magnitude  than  A  then  the  same  must  be  true  of  all  intermediate 
results  A/k ,  B/m ,  C/m  and  D/k  ,  so  the  final  results  are  right. 

function  Quotient  ( A ,  S,  C,  D)  :  ...  to  get  (A/B) -r  (C/D'i  ilt 

return  Quotient  l~  Product < A,  B,  D,  C)  ; 
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logical  function  Equal (A,  B ,  C,  D)i  ...  does  A/B  =  C/P  ? 
if  (B=0  and  0=0)  then 

{  if  <4=0  or  C=0>  then  Equal  :=  FALSE 
else  Equal  l-  TRUE  3 

else  <.  if  (A—C  and  B=D)  then  Equal  :=  TRUE 

else  Equal  : =  FALSE  }  s 

return  Equal  ; 
end. 

This  procedure  Equal  depends  crucially  upon  the  presumption  that 
its  arguments  A/B  and  C/D  are  ilt.  Note  also  that  0/0  is  not 
equal  to  anything,  not  even  itself,  since  it's  “Not  a  Number." 

Addition  and  subtraction  are  complicated  procedures  because  thev 
have  to  cope  with  expressions  like  4*P+G*C  when  their  values  are 
somewhat  smaller  in  magnitude  than  A  even  though  the  individual 
products  are  not.  The  following  subprocedure  is  needed. 


Coping  with  the  Determinant  -  y *z  : 

The  evaluation  of  expressions  like  x*t  -  y *z  =  det(jr  f)  when 
x ,  y,  z ,  t  and  the  determinant  are  all  integers  somewhat  srnal  J  er 
than  A  in  magnitude,  even  though  and  y*z  are  both  rather 

bigger,  is  a  subtask  that  occurs  often  enough  to  deserve  separate 
attention.  Our  approach  is  inspired  by  Gaussian  Elimination 
except  that,  instead  of  seeking  a  biggest  pivot  in  order  to 
secure  numerical  stability,  it  finds  the  smallest  element  in  the 
array  (;  j)  and  reduces  some  other  element  to  half  that  sire. 

The  reduction  process  ends  either  when  x*t  and  y*z  differ  in 
sign,  or  when  they  are  both  smaller  than  A  ,  in  which  cases  the 
determinant  can  be  evaluated  safely. 


whi  1 1 
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n  ‘ =  integer  nearest  y/x  ; 
y  )=  y  -  x*n  ;  ...  =  remfy.  x> 

t  1=  t  -  z*n  ;  ...  =  (Det  +  y*z)/x 

...  now  (new  y  I  <_  |.v/2|  and 

...  I  new  t|  <  IPet/.vl  +  lr/21  . 


return  Det  :=  -  y*r  : 

end . 
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Addition  and  Subtraction: 

Like  the  -foregoing  -functions 
procedures  act  upon  two  pairs 


Product  and 
o-f  integers 


Quotient , 
that  will 


the 

pass 


•foil  owi  ng 
unchanged 


through  the  -function  lit,  and  the  results  are 
likewise  provided  their  magnitudes  are  somewhat 


pairs  that  will 
less  than  A  . 


do 


function  Sum(A ,  B,  C ,  B)i  ...  to  get 
return  Sum  :=  DiffM,  S,  -C,  D)  ; 
end. 


( A/B )  +  (C/D)  ilt 


i  It 


function  Biff (A,  B ,  C,  D) :  ...  to  get  (A/B)  -  (C/D) 

G  max { 1 ,  gcd  <B , D) J  ;  b  1*  B/G  ;  d  t=  D/G  ; 

Now  we  seek  (A*d-b*C) / (G*b*d)  ilt,  but  first  we 
must  cancel  any  common  factor  g  hiding  in  G  : 
rem(A,G.)  ;  c  I  =  rem(C,G)  t  g  ‘.=  gcd(G,  a*d-b*c )  ; 
Note  | a*d-b*c I  <  | d*G/2 I + I b*G/2 I  <  A  . 

H  (G/g)*b*d  ;  ...  the  desired  denominator. 

The  numerator  wi 1 1  be  M  —  (A*d-b*C) /g  ... 
a  :=  rem(a,g.)  ;  c  C=  rem(c,gJ  ; 

H  :=  (a*d-b*c) /g  +  Det (  ( A-a)/g ,  b ,  (C-c)/g,  d  )  ; 

...  Note  how  |a*d-b*c|  <  A  as  before, 
return  Diff  (W,  N)  ; 
end. 


J  Are  they  worth  the  bother? 

'  It  seems  at  first  unlikely  that  a  calculation  of 

I  H/N  :=»  A/B  ±  C/D  *  (/!*£)  +  B*C)  /  (B*D)  ilt 


would  start  with  integers  A ,  B,  C,  D  not  much  smaller  than  A 
and  end  with  integers  M,  N  no  bigger  than  A  .  But.  having 
programmed  the  foregoing  procedures  into  various  prograrnmabl  e 
calculators  including  an  HP-97  and  an  HP-71 B,  I  have  seen 
these  unlikely  events  occur  about  as  often  as  not.  Perhaps  this 
is  merely  evidence  that  I  have  been  computing  some  things  the  hard 
way  instead  of  the  easy  way,  rather  than  evidence  that  anyone 
else  will  use  the  programs  every  day. 


These  programs  are  the  simplest  I  know  that  exemplify  a  property 
more  often  found  among  numerical  programs  than  others:  their 
simple  and  natural  specifications  belie  complicated  and  unnatural 
implementations.  It  may  seem  natural  to  demand  that,  if  the  dat^ 
given  a  program  and  the  output  desired  from  it  can  both  be 
represented  exactly  within  the  convenient  range  of  a  computer's 
capabilities,  the  output  actual ly  delivered  should  be  correct. 

But  that  demand  implies  that  the  program  will  find  a  path  from  the 
data  to  the  output  without  first  tr ansgr essi ng  the  computer's 
limitations  despite  that  the  path  begins  and  ends  only  a  step  or 
two  away  from  the  edge.  Such  a  path  need  not  be  obvious. 


Programs: 

Programs  for  the  HP-67/97  and  HP-71 B  have  been  appended  to 
these  notes.  The  program  for  the  HP-67/97  requires  very  little 
change  to  run  on  the  HP-4 1C  or  HP-15C.  Although  the  HP-7  IB 
program  is  written  in  a  kind  of  BASIC  that  looks  as  if  it  would 
run  on  diverse  other  machines,  the  program  exploits  the  HF-71B  s 
conformity  to  IEEE  p854  in  two  ways.  First,  its  rern  operator 
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(called  RED  on  the  HP-71B)  is  built-in  and  allows  the  program 
to  handle  integer  inputs  as  big  as  A  =  100  00000  00000  .  Second, 
the  Inexact  signal  accessible  through  FLAG(INX,  ...)  permits 
the  program  to  try  obvious  algorithms  first  and  then,  only  if  it 
encounters  roundoff,  resort  to  slower  ones.  Chained  sequences  of 
rational  operations  can  be  attempted  in  confidence  because  their 
results  will  assuredly  be  correct  unless  Inexact  is  signaled. 
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HP-67/97  prograa  to  perfora  RATIONAL  ARITHMETIC  on  pairs  of  integers  in  Lowest  Tern 

Usage:  The  stack  holds  four  integers  I,  Y,  7,  T  construed  as  two  rational  nuobers  Y/X  and  T/7  , 

both  presuaed  to  be  in  lowest  teras  (ilt).  If  not,  pressing  IE1  will  reduce  Y/X  to  lowest 
tern  while  leaving  T/Z  unchanged.  The  four  rational  operations  are  perforaed  by  pressing  one  of 
the  keys  [A],  IB],  (Cl,  101  to  invoke  reliable  progran,  or  Cal,  lb]  Cc)  Cdl  to  invoke  obvious 
program.  The  reliable  prograos  accept  integers  as  large  as  1,999,999.999  and  deliver  exactly 
correct  results  up  to  8,000,000,000  .  Specifically,  the  prograas  ... 

Add:  Press  CA1  or  Cal  to  put  Y/X  :=  (T/Z)  ♦  (Y/X)  ilt,  leaving  T/Z  unchanged. 

Subtract:  Press  CB1  or  lb]  to  put  Y/X  1*  (T/Z)  -  (Y/X)  ilt,  leaving  T/Z  unchanged. 

Multiply:  Press  CC1  or  tel  to  put  Y/X  Is  (T/Z)  *  (Y/X)  ilt,  leaving  T/Z  unchanged. 

Divide:  Press  ID!  or  Cd]  to  put  Y/X  1*  (T/Z)  f  (Y/X)  ilt,  leaving  T/Z  unchanged. 

Reduce:  Press  CEJ  to  put  Y/X  (Y/X)  ilt,  leaving  T/Z  unchanged. 

SCO:  Press  Cel  to  put  X  Greatest  Coaaon  Divisor  of  X  and  Y  . 

REN:  Press  CSSB1  C81  to  put  X  I*  Y  -  nX  and  n  :=  Integer  nearest  Y/X  into  reg.  8  . 

The  prograas  use  registers  0  to  8  and  I,  and  labels  2  to  8  too. 

Prograa:  *18L  A  CHS  *181  B  SSB  7  XJY  Rf  STO  4  6SB  e  X=0?  EEX  STO  5  STOfO  STOff  RCL  I  i<>Y 
6SB  8  RCL  «  x  STO  6  RCL  3  Rf  6SB  8  RCL  0  STO  7  x  RCL  6  -  6SB  e  ST0f5  RCL  I  X<>Y 

SSB  8  RCL  «  STOxO  x  STO  1  RCL  8  STO  6  RCL  3  Rf  SSB  8  RCL  7  x  RCL  t  -  W  f  STO  I 

RCL  5  STOxO  RCL  8  STO  5  »LBL  5  RCL  6  RCL  4  x  ENTt  ENTt  RCL  5  RCL  7  x  x  EEX  1  0 

X>Y?  STO  ♦  RCL  6  ABS  RCL  4  ABS  X<Y?  STO  3  LASTX  RCL  6  STO  4  W  STO  6  *L8L  3  RCL  5 
ABS  RCL  7  ABS  X<Y?  STO  3  LASTX  RCL  5  STO  7  X7Y  STO  S  *LBL3  RCL  7  ABS  RCL  4  ABS 
X<Y?  6T0  3  RCL  6  RCL  7  SSB  8  STO  6  RCL  4  RCL  8  x  STO-5  STO  5  *LBL  3  RCL  5  RCL  4 

6SB  8  STO  3  RCL  7  RCL  8  x  STO-6  STO  5  *LBL  4  LASTX  Rt  -  STO*!  6T0  6 

♦LSI  e  3  CHS  STO  I  Rt  1=0?  STO  2  SSB  8  6T0  i  (  juaps  back  three  steps  to  X=<T  ) 

•LBL  8  STO  3  X*Y  ENTt  ENTt  RCL  8  t  DSP  0  RND  STO  8  Rf  x  -  RTN  HBL  2  W  ABS  RTN 

♦LBL  7  STO  0  Rt  STO  I  Rt  STO  2  Rt  STO  3  RTN 

•LBL  E  SSB  7  Rt  6SB  e  X=0?  STO  6  STOfO  STOf ! 

•LBL  6  RCL  3  RCL  2  RCL  I  RCL  0  X>0?  RTN  CHS  I*Y  CHS  W  RTN 

•LBL  D  X*Y  »LBL  C  6SB  7  STO  4  SSB  e  X40?  STOfO  140?  ST0f4  RCL  I  RCL  2  SSB  e 

X40?  STOtl  RCL  2  X*Y  X40?  i  STOxO  RCL  4  STOxI  STO  6 

♦LBL  a  CHS  »LBL  b  SSB  7  x  X*Y  Rt  STOxO  x  -  STO  1  SS B  6  STO  E 

•LBL  d  XtY  »LBL  c  6SB  7  STOxl  Rt  STOxO  SSB  6  STO  E 
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10  !  Listing  of  HP-71B  program  to  perform  RATIONAL  ARITHMETIC 
20  !  upon  pairs  of  integers  in  Lowest  Terms  conveyed  as 
30  "Complex  Variables"  to  represent  R  =  M/N  as  <M,N)  . 

40  !  The  "Complex"  functions  herein  are  ... 

50  !  f nA  <R , S)  =  R+S  fnS<R,S)  =  R-S 

60  !  fnM(R,S)  =  R*S  fnD(R,S>  =  R/S 

70  !  fnl(R)  *  R  in  lowest  terms  (ilt) 

SO  !  Supporting  Real  functions  include  ... 

90  !  fnD0(R,S>  =  det(R,S)  =  Irnpt (Con j (R) *S) 

100  !  fnS(I,J)  =  Greatest  Common  Divisor  of  I  and  J  . 

110  !  RED  < I , J )  =  rem(I.J)  =  I  rem  J  as  in  IEEE  st ' d  p854 

120  !  RUN  to  sense  FLAG (INX)  and  reset  it  to  0  ;  if  that 
130  !  changes  then  a  result  has  been  compromised  by  roundoff. 
140  COMPLEX  R , S ,  R1,S1,  R2,S2,  R3,S3,  R4,S4,  R5,S5,  R6 
150  !  **■»***■»*****■»** 

160  DEF  FNG( 10, JO)  '  ...  =  GCD <10, JO) 

170  IF  J0=0  THEN  190 

ISO  o0=J0  @  J0=RED  <10,  JO)  <2  I0=o0  ®  IF  J0#0  THEN  180 
190  FNG= ABS  < 10)  <k  END  DEF 
200  !  ****•**•*•****#•*■*••* 

210  DEF  FNI(R6)  !  ...  =  R6  IN  LOWEST  TERMS 

220  FNI=R6/MAX (1 , FNG (REPT <R6>  , IMPT(R6) ) ) *SGN (CLASS ( I MPT (R6> )  ) 

230  END  DEF 

240  !  **•#-*■**■*■**■**■***-* 

250  DEF  FND0(R5,S5)  !  ...  =  det(R5,S5>  =  Irnpt (Con j (R5) *S5) 

260  ol =REPT (R5)  @  o2=IMPT  <R5)  @  o3=REPT(35)  &  o4=IMPT(S5) 

270  oO=FLAG ( I NX , 0 )  @  o5=SGN (ol *o4> *SGN <o2*o3>  @  oO=FLAG ( I NX , oO ) 
280  IF  o0=0  OR  o5#l  THEN  350 

290  IF  ABS ( o3 ) >ABS ( o2 )  THEN  o5*o3  @  o3=o2  @  o2=o5 

300  IF  ABS (a  1 ) >ABS  <o4)  THEN  o5=ol  @  ol*o4  @  o4=o5 

310  IF  ABS  < o 1 ) < = ABS  < o3 )  THEN  330 

320  a5=ol  @  a  1  =-o3  @  a3=o5  <3  o5=o2  @  o2=-o4  &  o4=o5 

330  o5=RED (o2,ol )  @  o0= (o2-o5) /ol  @  o2=o5  @  o4=o4-o3*o0 
340  GOTO  270 

350  FND0=ol*o4-o2*a3  @  END  DEF 
360  !  **•»**■»**#***#** 

370  DEF  FNM(R4,S4>  !  ...  =  R4-»S4  in  lowest  terms 

380  ol  =REPT  (R4)  @  o2=»IMPT(R4)  @  o3=REPT(S4)  @  o4=IMFT(S4) 

390  o5=MAX ( 1 , FNG ( a  1 , o4 ) )  @  oO=MAX ( 1  , FNG <a2 , o3) ) 

400  FNM=( (ol /o5) *(o3/o0) ,  (o2/o0) * (o4/o5) )  &  END  DEF 

410  !  ###**#*■**#**■*** 

420  DEF  FND <R3, S3) =FNM (R3, ( IMPT (S3) , REPT (S3) ) )  !  ...  =  R3/S3  lit 

430  !  ■»***■»****■»****■» 

440  DEF  FN S(R2,S2)  !  ...  =  R2-S2  in  lowest  terms 

450  ol  =REPT  (R2)  @  o2=IMF'T(R2)  0  o3=REPT(S2)  @  o4=IMFT(S2) 

460  oO=FLAG  ( INX  ,0)  @  o5=o  1  *o4-o2*o3  @  o6=o2*o4  <3  oO=FLAG  ( I  NX  ,  oO ) 

470  IF  o0=0  THEN  FNS=FNI  <  (o5,o6)  )  <2  GOTO  530 

480  o9=MAX ( 1 ,FNG(o2,o4) )  @  o2=o2/o9  &  o4=o4/o9 

490  o6=RED (o 1 , o9)  @  o7=RED (o3 . o9)  @  o5=FNG ( o9 . o6*o4-o2*o7) 

500  o9= (o9/o5) *o2*o4  &  o6=RED (06, o5)  @  o7=RED (o7 , o5) 

510  a8=  (a6*o4-o2*o7) /o5  @  o  1  =  (o  1 -06 ) /o5  (1  o3=  <o3-o7) /o5 

520  FNS=<FND0( (ol ,o2) , (o3,o4) )+o8,  o9) 

530  END  DEF 

540  '  ***************** 

550  DEF  FNA ( R 1 , S 1 ) =FNS ( R 1  ,  ( -REPT (31  )  , I MPT (SI )  >  >  '  ...  =  Rl+Sl  lit 
560  !  ***■**•*•»■»**»**** 

570  IF  FLAG ( INX , 0) =0  THEN  DISP  "Exact"  ELSE  DISF  "Inevact" 
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Calculating  M/N  I*  A/B  +  C/D  in  lowest  terms,  given 
integers  Af  5,  C  and  D,  is  a  task  taught  in  Elementary 
schools;  and  it  is  an  easy  exercise  in  Computer  Programming 
too  provided  the  given  integers  must  be  less  than  half  as  wide 
as  the  widest  integers  that  can  be  handled  conveniently  by  the 
computer's  hardware  or  by  its  programming  language.  But  that 
program  becomes  much  more  complicated  (and  slower)  if  it  is 
naively  expected  to  perform  correctly  whenever  all  six  of  our 
integers  A,  B ,  C,  D,  M  and  N  are  allowed  to  grow  almost  as 
wide  as  those  widest  convenient  integers.  This  simple  tas) 
illustrates  why  the  art  of  programming  entails  sometimes  a 
delicate  balance  between,  on  the  one  hand,  the  simplicity 
and  aesthetic  appeal  of  the  specifications  and,  on  the  other 
hand,  the  complexity  and  efficiency  of  the  implementation. 
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